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ABSTRACT 

We numerically investigate the dynamics of a 2D non-magnetised protoplanetary disc surrounded by an inflow coming from an 
external envelope. We find that the accretion shock between the disc and the inflow is unstable, leading to the generation of large- 
amplitude spiral density waves. These spiral waves propagate over long distances, down to radii at least ten times smaller than the 
accretion shock radius. We measure spiral-driven outward angular momentum transport with 10“'* < a < 10“^ for an inflow accretion 
rate > 10“* Mq ■ yr“*. We conclude that the interaction of the disc with its envelope leads to long-lived spiral density waves and 
radial angular momentum transport with rates that cannot be neglected in young non-magnetised protostellar discs. 
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1. Introduction 


Accretion is an essential phenomenon in astrophysics as it is the 
way through which gravitational energy is transformed into heat 
and therefore observable radiation. However, the physical pro¬ 
cess responsible for accretion is still debated. Although the mag- 
netorotational instability (MRI, |Balbus & Hawle^ 1991| l pro¬ 
vides an efficient accretion mechanism, its applicability to pro¬ 
toplanetary discs is questionable since these objects are very 
weakly ionised and might quench the MRI through various non¬ 
ideal magnetohydrodynamical (MHD) processes (ITumer et al.l 
2 ^ . 


The usual approach to protoplanetary disc dynamical mod¬ 
elling is to assume that these objects are isolated in space. In 
this context, most of the models rely on processes such as MRI- 
driven turbulence, self-gravity, or winds to drive accretion. How¬ 
ever, these discs are not isolated systems, and the question is 
whether surrounding material, albeit much less dense than the 
disc, could perturb it sufficiently to drive accretion. This possi¬ 
bility was proposed by |Padoan et al.| ( 2005[ l to explain accretion 
rates that scale like the central protostar mass squared. This sce¬ 
nario was later refined by Throo p & Bally| ( |2008| l, |Klessen &] 
Hennebelle (|2010[), and Padoan et al.|(2014)l and interpreted as a 


consequence of Bondi-Hoyle accretion onto the protostar. 


In this letter, we explore the importance of external accretion 
using a very simplified model. We consider a 2D hydrodynamic 
Keplerian disc onto which falls gas coming from a surrounding 
envelope. The disc is inviscid and hydrodynamically stable so 
that no accretion can occur without this external inflow. We first 
present in detail the model and the numerical method used to 
solve the equations of motion. We then explore some of the re¬ 
sults obtained using this model and finally discuss the limitations 
and implications of our findings. 


This work shares some similarities with that of IBae et al.l 
P015| l: both consider accretion discs subject to an external in¬ 
flow using 2D hydrodynamical models. However, we focus here 
on the transport of angular momentum far from the perturbed 
region using extended disc models, whereas |Bae et al.| ( |20T5] l 
concentrated on the local response of the disc in the presence of 
viscous angular momentum transport. We discuss the differences 
between the outcome of these models in Sect. 4. 


2. Model 


2.1. Physical model 


We study the dynamics of an accretion disc close to Keplerian 
rotation. We consider a 2D frame (R, 0) and integrate the equa¬ 
tions of motion in the vertical direction. In this initial study, we 
neglect the role played by magnetic fields and self-gravity to 
concentrate on a purely hydrodynamical model. The equations 
of motion read 


— + V . (Ed) = 0, 
ot 

dv 1 

— -i-d*Vd = - —VP- 


GMo 


( 1 ) 

( 2 ) 


where v is the flow velocity, E the surface density, P the verti¬ 
cally averaged pressure, G the gravitational constant, and Mq the 
mass of the central object, which is assumed to be a solar-mass 
star. 

Instead of the usual minimum mass solar nebula (MMSN) 
density profile, we use a power-law surface density profile with 
an exponential outer cut-off 


E = 





(3) 
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where we choose Sq = 1700 g.cm and Rc = 40a.u., which 
correst tonds to typical proto planetary discs observed in Ophi- 


uchus ( Andrews et al.|2069 i with a mass of 4.8 X 10 ^ Mq. The 


temperature profile is chosen to be locally isothermal so that 
Cj/vk = s = 0.1 or 0.05, where Cj is the sound speed and 
vk = ^IGMqIR is the Keplerian velocity. The equation of state 
is then simply P - Assuming a vertically isothermal hydro¬ 
static profile, this model leads to an accretion disc with a con¬ 
stant aspect ratio HjR - e, where H is the typical vertical disc 
scale height. 


2.2. Boundary conditions 

The outer radial boundary condition is located at /?out = 400 a.u. 
To mimic matter falling onto the disc, we inject material at the 
outer radial boundary condition that then falls onto the disc. The 
material is injected over an azimuthal extent 0 < 0 < 0inj with 
a sonic radial velocity. = -Cj. We vary three parameters 
for the injected material: its specific angular momentum Xout = 
ve(f?out)^out, its accretion rate = -f d0SoutVr„„„ Siout being 
deduced from the desired accretion rate, and the azimuthal width 
of the inflow di„j. 

The inner radial boundary is located at Rin - 1 a.u. The inner 
boundary condition is forced to the initial density and velocity 
profile. To avoid spurious reflexion of spiral density waves, we 
add a damping zone in the region 1 a.u < /? < 3 a.u. In this 
region, we exponentially relax non-axisymmetric velocity fluc¬ 
tuations on a fixed timescale set to eight local orbital periods. To 
avoid mass accumulation in the damping zone, we also relax the 
density to the initial density profile on the same timescale. 


2.3. Numericai aigorithm 


We use Pluto 4.0 ( [Mignone et al.|2007) l to solve the equations of 
motion on a cylindrical grid, using log-spaced grid cells in the 
radial direction. Parabolic reconstruction is used in each cell to 
cope with the non-uniform radial grid, and a third-order Runge- 
Kutta algorithm is used to evolve the system in time. We use 
the HLL Riemann solvetQ at cell boundaries combined with or¬ 
bital advection ( |Mignone et al.|20T2 1 to allow for long integra¬ 
tion times. We have tested that starting with sonic white noise 
perturbation on v and without any inflow, the disc was hydro- 
dynamically stable with a ~ 10“^ after a few local orbits, as 
expected with our density profile Q. The resolution of each run 
was (NR,Ng) = (512,512)fore = 0.1 oi:(Nii,Ng) = a024,1024) 
for s = 0.05, which corresponds to a locally isotropic resolution 
of eight points per H. We checked that doubling the resolution of 
run S7-1 did not quantitatively affect our results. We are there¬ 
fore confident that our simulations are numerically converged. 


2.4. Units, diagnostics, and notations 

Length units are given in a.u., surface densities in g.cm^^, time 
units in orbital periods at 100 a.u. - T loo (or equivalently in units 
of 10^ yr), and accretion rates in Mo.yr^*. 

In the following, several diagnostics are used to measure the 
behaviour of the disc coupled to the inflow. We first introduce 
two averages, an azimuthal average ~ and a temporal average 

* We also ran simulations with the HLLC Riemann solver. In this case, 
however, the noise level of the background spirals in the absence of in¬ 

flow is O' ~ 10“"* because there is almost no numerical dissipation. We 
therefore only present results obtained with HLL for which the back¬ 
ground noise level is much lower. 
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Run 

Mint 

-^out 

a(Rd/2) 

a(Rd/10) 

Rd 

S8-1 

10-** 

1 

1.6 

X 

10-3 

5.3 

X 

10-3 

255 

S7-1 

10-^ 

1 

4.5 

X 

10-3 

1.4 

X 

10-4 

199 

S6-1 

10-® 

1 

1.0 

X 

10-3 

3.5 

X 

10-4 

167 

A7-1 

10-^ 

1 

3.5 

X 

10-3 

3.8 

X 

10-4 

184 

S8-1-thin 

10-** 

1 

3.8 

X 

T(r^ 

1.3 

X 

10-3 

225 

S7-l-thin 

10-^ 

1 

8.3 

X 

10-4 

3.4 

X 

10-3 

164 

S6-1-thin 

10-® 

1 

1.5 

X 

10-3 

1.1 

X 

10-4 

135 

S9-0 

10*^ 

0 

8.1 

X 


2.8 

X 

10-3 

264 

S8-0 

10-** 

0 

2.3 

X 

10-3 

1.0 

X 

10-4 

177 

S7-0 

10-^ 

0 

4.2 

X 

10-3 

3.0 

X 

10-4 

91 

A7-0 

10-’ 

0 

4.5 

X 

10-3 

5.0 

X 

10-4 

90 


Table 1. List of simulations. Mm is measured in M^.yr ', and Lout in 
units of specific angular momentum of a Keplerian disc at R = 100: 
Lioo- 


(■). We start the temporal average from t - 10 Tioo to allow the 
system to relax from the initial condition (this corresponds to 
10^ orbits at the inner boundary). 

These averaging procedures allow us to define fl uctuating 
quantities 6v — v — 'Zvf'Z. The first diagnostic is the 


Shakura 


& Sunyaev (1973i a parameter a = 'ZdvRdvgIP. We also use 


the disc radius R^ defined as the location where the disc rotation 
velocity drops below 90% of the Keplerian velocity vk. This def¬ 
inition might look rather arbitrary, but it nevertheless leads to a 
well-defined disc radius since the rotation velocity drops very 
rapidly with radius at the transition region between the disc and 
the inflow. Each simulation is integrated for 3x10^ orbits at the 
inner boundary, which corresponds to 30 T mo. 

Our runs are labelled "XY-Z" where X=S or A denotes sym¬ 
metric (0inj = 2n) or asymmetric (0inj = 7r/2) inflows, Y is the 
accretion rate at Rout and Z is the specific angular momentum 
Lout- Our simulations all have e = 0.1, except for runs labelled 
"thin", which assume s = 0.05. 


3. Results 

We first concentrate on symmetric inflows. We then explore the 
impact of asymmetric inflows and of the disc thickness in a sec¬ 
ond part. 

3.1. Symmetric accretion 

3.1.1. Accretion shock 

The initial evolution of the symmetric accretion run S7-1 is pre¬ 
sented in Fig. [T] The falling material initially forms an accretion 
shock that propagates inward. At f ^ 2rioo, the shock stalls at 
R = Rd- The stationary shock then develops an instability, which 
appears as a short-wavelength non-axisymmetric ondulation of 
the shock front (f = 2.8 Tioo)- This instability then quickly satu¬ 
rates and generates strong spiral waves that propagate inward on 
the long term (f > 3.57’mo). 

The origin of this instability is tightly linked to the structure 
of the shock. Spirals are continuously produced from R < Rd, 
and the instability survives for the entire duration of the simula¬ 
tion. The instability mechanism becomes clear when the average 
specific angular momentum profile (X) and the average potential 
vorticity {^) s ((2Q -H RQ')/S) are plotted as a function of ra¬ 
dius (Fig. [^. We And that at the location of the stationary shock 
a strong velocity gradient develops, which leads to an extremely 
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Fig. 1. Density maps as a function of time for run S7-1. From left to right: t = 2rioo, t = 2.8 Tigo, t = 6 Fioo- 



Fig. 2. Averaged potential vorticity {() (plain blue) and specific angular 
momentum X. (red dashed) of run S7-1. Note the sharp negative angular 
momentum gradient associated with a strong potential vorticity peak at 
R ^ 200. 



Fig. 3. Profiles of a as a function of radius, measured from Rg. The 
spiral stress increases with increasing Mi„f and with decreasing initial 
angular momentum Xout- 


peaked f. This sharp structure appears because the specific an¬ 
gular momentum of the falling material is lower than that of the 
disc material at R = 200. A potential vorticity layer is therefore 
unavoidable if matter keeps falling onto the disc. 

As is well known, such a hydrodynamical structure natu¬ 
rally leads to a Kelvin-Helmholtz instability (KHl, also known as 
Rossby wave instability, or RWl, in the astrophysical contexj^. 
We note, however, that the shear rate is locally so strong that it 
also locally violates the Rayleigh criterion: in the shock region, 
the specific angular momentum is decreasing outward (Fig. |^. 
As stated above, this configuration is unavoidable given that the 
inflow has a lower angular momentum than the disc. The insta¬ 
bility driving the spiral is therefore a mixture of KHl (RWl) and 
Rayleigh centrifugal instability. It has a growth rate close to the 
orbital frequency at Rd and is sustained, for reasons stated above, 
during the entire duration of the simulation. 


3.1.2. Angular momentum transport 

Spiral waves generate a positive Reynolds stress through the 
entire disc, which we quantify as an a parameter (Fig. |^. We 
find that these spirals produce a very strong transport close to 
R ^ Rd with a ^ 1. They propagate inward from Rd and dissi¬ 
pate through shocks, reaching a > 10 “^ at Rd/2 and a > 10 
at Rd/10. A generation region (0.6Rd < R < Rd) can be defined 


^ In our case, the RWl is not due to a density bump but to a strong and 
localised shear. Nevertheless, the general RWl criterion, having a 


maximum (Lovelace et al. 1 


is satisfied. 


where the inflow mixes with the disc, spiral modes are excited, 
and a values are similar for all of our simulations. Likewise, 
there is a propagation region (R < Q.6Rd) where the disc den¬ 
sity structure is unaffected by the inflow except for propagating 
spiral waves. 

We see that a in the generation region is not affected by Mmf, 
but a in the propagation region (R < Q.6Rd) clearly is. More¬ 
over, for our strongest Mjnf (run S6-1 for instance), we observe 
a bump of transport for R ^ Q.6Rd, directly at the transition be¬ 
tween these two regions. This is because the mass of the gener¬ 
ation region has increased significantly due to the inflow while 
that of the propagation region has not, forming a jump in surface 
density at the interface. As a result, the generation region has a 
higher inertia and excites stronger waves at the interface with the 
propagation region. We note that smaller angular momentum in 
the inflow also leads to somewhat stronger spirals. These find¬ 
ings can be summarised by representing a at R - Rd/10 as a 
function of (Fig. [^. 

From these results, we can estimate that in the symmetric 
inflow case 

a{Rd/lQ). (4) 

yiO^^Mo -yr-'/ 

For inflows with significant rotation (-Cout = 1) we have ao - 
1.3 X 10“^ and y - 0.4, whereas for inflows without initial ro¬ 
tation we obtain ao = 3.0 x lO^'* and y - 0.5. These values 
constitute a lower bound for the average alpha found in the disc. 
Typical values of a at R = Rd/2 are 10 to 100 times higher than 
these estimates (Table [^. 
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Fig. 4. a measured = Rd/W as a function of Minf. Colours represent 
different Xom or H/R, diamonds correspond to symmetric inflow, and 
stars correspond to asymmetric accretion. The dashed line are fits given 
hy Eq. l|^. 


As is well known, standard accretion disc theory allows de¬ 
riving a theoretical mass accretion rate from the profile of a 
computed above ( [Balbus & Papaloizou| ( |1999| ), Eq. 28). We have 
checked that the measured accretion rates are consistent with this 
theoretical expectation, provided that M is computed on long 
enough averages (typically more than ten local orbital periods). 
This shows that standard accretion disc theory also applies to 
spiral-driven transport of angular momentum. This result does 


not imply, however, that “turbulent” heating follows S 

hakura 

& Sunyaev (1973i or-disc model (Balbus & Papaloizou 

19991 

This point cannot be tested in our model since we used a 

ocally 


isothermal equation of state. 


3.2. Parameter survey 


Since spiral waves are primarily sound waves, a dependence on 
the sound speed and therefore on H/R = e is expected. To test 
this effect, we reduce H/R by a factor 2 with s - 0.05. These 
simulations, labelled “thin”, show that the amplitude of the spi¬ 
rals is significantly reduced when H/R is reduced. We find a typ¬ 
ical attenuation of a factor 3 between e = 0.1 and e = 0.05. This 
is similar to the case of tidally induced spiral waves ( |Savonije1 
et al. 1994| l, despite an excitation mechanism different from ours. 
This similitude indicates that the wave propagation is most likely 
mainly responsible for this dependence. 

All the results presented above were obtained in a situation 
where gas is falling onto the disc in an axisymmetric way. To 
test the effect of this assumed symmetry on the results, we con¬ 
sider the case of accretion coming from a narrower stream (runs 
labelled A*). We find that the disc shape is largely unaffected 
for the inflow rates we used and produces spiral waves similar 
to the one found in the symmetric case (Fig.[T]). In particular, we 
do not find any non-axisymmetric feature close to the location 
where the stream hits the disc. In addition, the azimuthally aver¬ 
aged profiles of {^) and {X.} are very similar to the one found in 
the symmetric accretion case. The average flow is still unstable 
to the same instability as in the symmetric case. It is therefore 
expected to find that asymmetric accretion leads to a(R) profiles 
that are qualitatively comparable to those found in the symmetric 
case. 

By comparing a(r) obtained in the asymmetric case to the 
symmetric case, we find that asymmetric accretion in general 
drives a slightly more efficient angular momentum transport. 
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probably thanks to the excitation of low m modes of the insta¬ 
bility by the stream that propagates more easily. At R - Rd/^0, 
we obtain an a about two to four times larger than the one ob¬ 
tained from symmetric accretion simulations (Fig.|^. This sug¬ 
gests that the symmetric accretion scenario given by Eq.[^ con¬ 
stitutes a lower bound to spiral-driven angular momentum trans¬ 
port. 


4. Discussion 


We explored the impact of an external inflow on a protoplane¬ 
tary disc using 2D hydrodynamical simulations. We found that 
this interaction leads to the generation of strong spiral waves 
that propagate on long distances (typically down to radii smaller 
that Rd/lO). The resulting a at Rd/2 are larger than 10“^ and 
reach a few 10 at f?d/10 for inflow mass rates higher than 
Mo.yr"'. This stress varies strongly with Minf and H/R, but 
is not significantly affected by the geometry of the inflow. 

As shown in the introduction, this study is similar in nature 
to the work of |Bae et ak ( |2015| l. However, several important dif¬ 
ferences should be emphasised. First, our simulations only as¬ 
sumed a radial inflow and no vertical inflow. Second, the angu¬ 
lar momentum of the falling material is always smaller than that 
of the outer disc so that a radial shear layer is always present, a 
possibility only explored in one simulation by Bae et al. (2015 i. 
Because of these differences, we do not observe the formation of 
vortices, most probably because the instability at Rd involves a 
mixture of the RWI and of Rayleigh centrifugal instability thanks 
to the strong shear layer. We also obtained higher a values than 
did |Bae et ak] ( |2015| l in the generation region, probably due to 
the same shear layer. 

Finally, spiral waves such as we discussed here are now de¬ 
tectable with infrared continuum observations. Several authors 
have already reported direct observations of spiral patterns at the 
outer edge of evolved protostellar discs (e.g. |Muto et ar]|2012t 
|Benisty et al.|20l5] ). Whether an external inflow might be at the 
origin of these spiral pattern is an open question that we defer to 
a future publication. 
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